#pragma once

#include "vec4.h"
#include "box3.h"
#include "../arch/avx.h"

OPEN_JLIB_NS


struct Mat4 {

  f32v4 c0, c1, c2, c3;

  inline Mat4 () noexcept = default;
  inline Mat4 (f32 s00, f32 s10, f32 s20, f32 s30, f32 s01, f32 s11, f32 s21, f32 s31, f32 s02, f32 s12, f32 s22, f32 s32, f32 s03, f32 s13, f32 s23, f32 s33) noexcept
  : c0(f32v4{s00, s01, s02, s03})
  , c1(f32v4{s10, s11, s12, s13})
  , c2(f32v4{s20, s21, s22, s23})
  , c3(f32v4{s30, s31, s32, s33}) {
  }
  inline Mat4 (f32v4 c0, f32v4 c1, f32v4 c2, f32v4 c3) noexcept : c0(c0), c1(c1), c2(c2), c3(c3) {
  }
  inline Mat4 (Vec4 c0, Vec4 c1, Vec4 c2, Vec4 c3) noexcept : c0(c0.p), c1(c1.p), c2(c2.p), c3(c3.p) {
  }
  inline Mat4 (const f32* m) noexcept {
    asm volatile (R"(
    ld1 { %[dt0].4s }, [%[mem]], #16
    ld1 { %[dt1].4s }, [%[mem]], #16
    ld1 { %[dt2].4s }, [%[mem]], #16
    ld1 { %[dt3].4s }, [%[mem]])"
    : [dt0]"=x"(c0), [dt1]"=x"(c1), [dt2]"=x"(c2), [dt3]"=x"(c3), [mem]"+r"(m));
  }


  template<uwl Index>
  inline f32 c () const noexcept;


  inline Vec4 operator* (Vec4 v) const noexcept {
    f32v4 r;
    asm volatile (R"(
    fmul %[dst].4s, %[mt0].4s, %[vec].s[0]
    fmla %[dst].4s, %[mt1].4s, %[vec].s[1]
    fmla %[dst].4s, %[mt2].4s, %[vec].s[2]
    fmla %[dst].4s, %[mt3].4s, %[vec].s[3])"
    : [dst]"=&x"(r): [mt0]"x"(c0), [mt1]"x"(c1), [mt2]"x"(c2), [mt3]"x"(c3): [vec]"x"(v.p));
    return r;
  }
  inline Mat4 operator* (const Mat4& m) const noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    fmul %[dt0].4s, %[ma0].4s, %[mb0].s[0]
    fmul %[dt1].4s, %[ma0].4s, %[mb1].s[0]
    fmul %[dt2].4s, %[ma0].4s, %[mb2].s[0]
    fmul %[dt3].4s, %[ma0].4s, %[mb3].s[0]
    fmla %[dt0].4s, %[ma1].4s, %[mb0].s[1]
    fmla %[dt1].4s, %[ma1].4s, %[mb1].s[1]
    fmla %[dt2].4s, %[ma1].4s, %[mb2].s[1]
    fmla %[dt3].4s, %[ma1].4s, %[mb3].s[1]
    fmla %[dt0].4s, %[ma2].4s, %[mb0].s[2]
    fmla %[dt1].4s, %[ma2].4s, %[mb1].s[2]
    fmla %[dt2].4s, %[ma2].4s, %[mb2].s[2]
    fmla %[dt3].4s, %[ma2].4s, %[mb3].s[2]
    fmla %[dt0].4s, %[ma3].4s, %[mb0].s[3]
    fmla %[dt1].4s, %[ma3].4s, %[mb1].s[3]
    fmla %[dt2].4s, %[ma3].4s, %[mb2].s[3]
    fmla %[dt3].4s, %[ma3].4s, %[mb3].s[3])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [ma0]"x"(c0), [ma1]"x"(c1), [ma2]"x"(c2), [ma3]"x"(c3), [mb0]"x"(m.c0), [mb1]"x"(m.c1), [mb2]"x"(m.c2), [mb3]"x"(m.c3));
    return {r0, r1, r2, r3};
  }


  inline f32 operator* () const noexcept {
    f32 r, a0, b0, ab, t0, t1;
    asm volatile (R"(
    dup %[tp0].4s, %[mt0].s[0]                   // [s00, s00, s00, s00]
    ins %[tp0].s[3], %[mt0].s[1]                 // [s00, s00, s00, s01]
    ext %[tp1].16b, %[mt1].16b, %[mt1].16b, #4   // [s11, s12, s13, s10]
    ins %[tp1].s[3], %[mt1].s[3]                 // [s11, s12, s13, s13]
    fmul %[ta0].4s, %[tp0].4s, %[tp1].4s         // [a0, a1, a2, a4]
    ext %[tp0].16b, %[mt0].16b, %[mt0].16b       // [s01, s02, s03, s01]
    ins %[tp0].s[3], %[mt0].s[3]                 // [s01, s02, s03, s03]
    dup %[tp1].4s, %[mt1].s[0]                   // [s10, s10, s10, s10]
    ins %[tp1].s[3], %[mt1].s[1]                 // [s10, s10, s10, s11]
    fmls %[ta0].4s, %[tp0].4s, %[tp1].4s         // [a0, a1, a2, a4]
    dup %[tp0].4s, %[mt2].s[1]                   // [s21, s21, s21, s21]
    ins %[tp0].s[0], %[mt2].s[2]                 // [s22, s21, s21, s21]
    ins %[tp0].s[3], %[mt2].s[0]                 // [s22, s21, s21, s20]
    dup %[tp1].2s, %[mt3].s[3]                   // [s33, s33, -, -]
    ins %[tp1].s[2], %[mt3].s[2]                 // [s33, s33, s32, -]
    ins %[tp1].s[3], %[mt3].s[2]                 // [s33, s33, s32, s32]
    fmul %[tb0].4s, %[tp0].4s, %[tp1].4s         // [b0, b1, b2, b4]
    dup %[tp0].2s, %[mt2].s[3]                   // [s23, s23, -, -]
    ins %[tp0].s[2], %[mt2].s[2]                 // [s23, s23, s22, -]
    ins %[tp0].s[3], %[mt2].s[2]                 // [s23, s23, s22, s22]
    dup %[tp1].4s, %[mt3].s[1]                   // [s31, s31, s31, s31]
    ins %[tp1].s[0], %[mt3].s[2]                 // [s32, s31, s31, s31]
    ins %[tp1].s[3], %[mt3].s[0]                 // [s32, s31, s31, s30]
    fmls %[tb0].4s, %[tp0].4s, %[tp1].4s         // [b0, b1, b2, b4]
    ext %[tp0].16b, %[mt0].16b, %[mt2].16b, #4   // [s01, s02, s03, s20]
    ins %[tp0].s[2], %[mt2].s[0]                 // [s01, s02, s20, s20]
    ext %[tp1].16b, %[mt1].16b, %[mt3].16b, #8   // [s12, s13, s30, s31]
    ins %[tp1].s[2], %[mt3].s[3]                 // [s12, s13, s33, s31]
    fmul %[tab].4s, %[tp0].4s, %[tp1].4s         // [a3, a5, b3, b5]
    ext %[tp0].16b, %[mt0].16b, %[mt2].16b, #8   // [s02, s03, s20, s21]
    ins %[tp0].s[2], %[mt2].s[3]                 // [s02, s03, s23, s21]
    ext %[tp1].16b, %[mt1].16b, %[mt3].16b, #4   // [s11, s12, s13, s30]
    ins %[tp1].s[2], %[mt3].s[0]                 // [s11, s12, s30, s30]
    fmls %[tab].4s, %[tp0].4s, %[tp1].4s         // [a3, a5, b3, b5]
    ins %[tp1].d[0], %[tab].d[1]                 // [b3, b5, -, -]
    fmul %[tp1].2s, %[tab].2s, %[tp1].2s         // [a3b3, a5b5, -, -]
    eor %[tp0].8b, %[tp0].8b, %[tp0].8b          // [0, 0, -, -]
    zip1 %[tp1].2d, %[tp1].2d, %[tp0].2d         // [a3b3, a5b5, 0, 0]
    fmul %[tp0].4s, %[ta0].4s, %[tb0].4s         // [a0b0, a1b1, a2b2, a4b4]
    uzp1 %[det].4s, %[tp0].4s, %[tp1].4s         // [a0b0, a2b2, a3b3, 0]
    uzp2 %[tp1].4s, %[tp0].4s, %[tp1].4s         // [a1b1, a4b4, a5b5, 0]
    fneg %[tp1].2s, %[tp1].2s                    // [-a1b1, -a4b4, a5b5, 0]
    fadd %[det].4s, %[det].4s, %[tp1].4s         // [a0b0-a1b1, a2b2-a4b4, a3b3+a5b5, 0]
    faddp %[det].4s, %[det].4s, %[det].4s
    faddp %s[det], %[det].2s)"
    : [det]"=x"(r), [ta0]"=&x"(a0), [tb0]"=&x"(b0), [tab]"=&x"(ab), [tp0]"=&x"(t0), [tp1]"=&x"(t1)
    : [mt0]"x"(c0), [mt1]"x"(c1), [mt2]"x"(c2), [mt3]"x"(c3));
    return r;
  }


  inline Mat4 operator~ () const noexcept {
    f32v4 r0, r1, r2, r3, t0, t1;
    asm volatile (R"(
    zip1 %[tp0].2d, %[mt0].2d, %[mt2].2d         // [s00, s01, s20, s21]
    zip1 %[tp1].2d, %[mt1].2d, %[mt3].2d         // [s10, s11, s30, s31]
    trn1 %[dt0].4s, %[tp0].4s, %[tp1].4s         // [s00, s10, s20, s30]
    trn2 %[dt1].4s, %[tp0].4s, %[tp1].4s         // [s01, s11, s21, s31]
    zip2 %[tp0].2d, %[mt0].2d, %[mt2].2d         // [s02, s03, s22, s23]
    zip2 %[tp1].2d, %[mt1].2d, %[mt3].2d         // [s12, s13, s32, s33]
    trn1 %[dt2].4s, %[tp0].4s, %[tp1].4s         // [s02, s12, s22, s32]
    trn2 %[dt3].4s, %[tp0].4s, %[tp1].4s         // [s03, s13, s23, s33]
    )"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=x"(r2), [dt3]"=x"(r3), [tp0]"=&x"(t0), [tp1]"=&x"(t1)
    : [mt0]"x"(c0), [mt1]"x"(c1), [mt2]"x"(c2), [mt3]"x"(c3));
    return {r0, r1, r2, r3};
  }
  inline Mat4 operator- () const noexcept {
    f32v4 det, t0, t1, a0, b0, ab, r0, r1, r2, r3, s0, s1, s2, s3;
    asm volatile (R"(
    dup %[tp0].4s, %[mt0].s[0]                   // [s00, s00, s00, s00]
    ins %[tp0].s[3], %[mt0].s[1]                 // [s00, s00, s00, s01]
    ext %[tp1].16b, %[mt1].16b, %[mt1].16b, #4   // [s11, s12, s13, s10]
    ins %[tp1].s[3], %[mt1].s[3]                 // [s11, s12, s13, s13]
    fmul %[ta0].4s, %[tp0].4s, %[tp1].4s         // [a0, a1, a2, a4]
    ext %[tp0].16b, %[mt0].16b, %[mt0].16b       // [s01, s02, s03, s01]
    ins %[tp0].s[3], %[mt0].s[3]                 // [s01, s02, s03, s03]
    dup %[tp1].4s, %[mt1].s[0]                   // [s10, s10, s10, s10]
    ins %[tp1].s[3], %[mt1].s[1]                 // [s10, s10, s10, s11]
    fmls %[ta0].4s, %[tp0].4s, %[tp1].4s         // [a0, a1, a2, a4]
    dup %[tp0].4s, %[mt2].s[1]                   // [s21, s21, s21, s21]
    ins %[tp0].s[0], %[mt2].s[2]                 // [s22, s21, s21, s21]
    ins %[tp0].s[3], %[mt2].s[0]                 // [s22, s21, s21, s20]
    dup %[tp1].2s, %[mt3].s[3]                   // [s33, s33, -, -]
    ins %[tp1].s[2], %[mt3].s[2]                 // [s33, s33, s32, -]
    ins %[tp1].s[3], %[mt3].s[2]                 // [s33, s33, s32, s32]
    fmul %[tb0].4s, %[tp0].4s, %[tp1].4s         // [b0, b1, b2, b4]
    dup %[tp0].2s, %[mt2].s[3]                   // [s23, s23, -, -]
    ins %[tp0].s[2], %[mt2].s[2]                 // [s23, s23, s22, -]
    ins %[tp0].s[3], %[mt2].s[2]                 // [s23, s23, s22, s22]
    dup %[tp1].4s, %[mt3].s[1]                   // [s31, s31, s31, s31]
    ins %[tp1].s[0], %[mt3].s[2]                 // [s32, s31, s31, s31]
    ins %[tp1].s[3], %[mt3].s[0]                 // [s32, s31, s31, s30]
    fmls %[tb0].4s, %[tp0].4s, %[tp1].4s         // [b0, b1, b2, b4]
    ext %[tp0].16b, %[mt0].16b, %[mt2].16b, #4   // [s01, s02, s03, s20]
    ins %[tp0].s[2], %[mt2].s[0]                 // [s01, s02, s20, s20]
    ext %[tp1].16b, %[mt1].16b, %[mt3].16b, #8   // [s12, s13, s30, s31]
    ins %[tp1].s[2], %[mt3].s[3]                 // [s12, s13, s33, s31]
    fmul %[ta1].4s, %[tp0].4s, %[tp1].4s         // [a3, a5, b3, b5]
    ext %[tp0].16b, %[mt0].16b, %[mt2].16b, #8   // [s02, s03, s20, s21]
    ins %[tp0].s[2], %[mt2].s[3]                 // [s02, s03, s23, s21]
    ext %[tp1].16b, %[mt1].16b, %[mt3].16b, #4   // [s11, s12, s13, s30]
    ins %[tp1].s[2], %[mt3].s[0]                 // [s11, s12, s30, s30]
    fmls %[tab].4s, %[tp0].4s, %[tp1].4s         // [a3, a5, b3, b5]
    ins %[tp1].d[0], %[tab].d[1]                 // [b3, b5, -, -]
    fmul %[tp1].2s, %[tab].2s, %[tb1].2s         // [a3b3, a5b5, -, -]
    eor %[tp0].8b, %[tp0].8b, %[tp0].8b          // [0, 0, -, -]
    zip1 %[tp1].2d, %[tp1].2d, %[tp0].2d         // [a3b3, a5b5, 0, 0]
    fmul %[tp0].4s, %[ta0].4s, %[tb0].4s         // [a0b0, a1b1, a2b2, a4b4]
    uzp1 %[det].4s, %[tp0].4s, %[tp1].4s         // [a0b0, a2b2, a3b3, 0]
    uzp2 %[tp1].4s, %[tp0].4s, %[tp1].4s         // [a1b1, a4b4, a5b5, 0]
    fneg %[tp1].2s, %[tp1].2s                    // [-a1b1, -a4b4, a5b5, 0]
    fadd %[det].4s, %[det].4s, %[tp1].4s         // [a0b0-a1b1, a2b2-a4b4, a3b3+a5b5, 0]
    faddp %[det].4s, %[det].4s, %[det].4s
    faddp %[det].4s, %[det].4s, %[det].4s
    fneg %[dt0].4s, %[mt0].4s                    // [-s00, -s01, -s02, -s03]
    fneg %[dt2].4s, %[mt2].4s                    // [-s20, -s21, -s22, -s23]
    zip1 %[tp0].2d, %[dt0].2d, %[dt2].2d         // [-s00, -s01, -s20, -s21]
    zip1 %[tp1].2d, %[mt1].2d, %[mt3].2d         // [s10, s11, s30, s31]
    trn1 %[sm0].4s, %[tp1].4s, %[tp0].4s         // [s10, -s00, s30, -s20]
    trn2 %[sm1].4s, %[tp1].4s, %[tp0].4s         // [s11, -s01, s31, -s21]
    zip2 %[tp0].2d, %[dt0].2d, %[dt2].2d         // [-s02, -s03, -s22, -s23]
    zip2 %[tp1].2d, %[mt1].2d, %[mt3].2d         // [s12, s13, s32, s33]
    trn1 %[sm2].4s, %[tp1].4s, %[tp0].4s         // [s12, -s02, s32, -s22]
    trn2 %[sm3].4s, %[tp1].4s, %[tp0].4s         // [s13, -s03, s33, -s23]
    eor %[dt1].16b, %[dt1].16b, %[dt1].16b       // [0, 0, 0, 0]
    mov %[dt3].16b, %[dt1].16b                   // [0, 0, 0, 0]
    rev64 %[tab].2s, %[tab].2s                   // [a5, a3, b3, b5]
    zip1 %[tp0].2s, %[tb0].2s, %[tab].2s         // [b0, a5, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b0, b0, a5, a5]
    fmul %[dt0].4s, %[sm1].4s, %[tp0].4s
    fmls %[dt1].4s, %[sm0].4s, %[tp0].4s
    ext %[tp1].16b, %[ta0].16b, %[ta0].16b, #8   // [a2, a4, a0, a1]
    zip2 %[tp0].2s, %[tb0].2s, %[tp1].2s         // [b1, a4, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b1, b1, a4, a4]
    fmls %[dt0].4s, %[sm2].4s, %[tp0].4s
    fmul %[dt2].4s, %[sm0].4s, %[tp0].4s
    ext %[tb0].16b, %[tb0].16b, %[tb0].16b, #4   // [b1, b2, b4, b0]
    zip2 %[tp0].2s, %[tb0].2s, %[tab].2s         // [b2, a3, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b2, b2, a3, a3]
    fmla %[dt0].4s, %[sm3].4s, %[tp0].4s
    fmls %[dt3].4s, %[sm0].4s, %[tp0].4s
    ext %[tab].16b, %[tab].16b, %[tab].16b, #8   // [b3, b5, a5, a3]
    zip1 %[tp0].2s, %[tab].2s, %[tp1].2s         // [b3, a2, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b3, b3, a2, a2]
    fmla %[dt1].4s, %[sm2].4s, %[tp0].4s
    fmls %[dt2].4s, %[sm1].4s, %[tp0].4s
    ext %[tb0].16b, %[tb0].16b, %[tb0].16b, #4   // [b2, b4, b0, b1]
    zip2 %[tp0].2s, %[tb0].2s, %[ta0].2s         // [b4, a1, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b4, b4, a1, a1]
    fmls %[dt1].4s, %[sm3].4s, %[tp0].4s
    fmla %[dt3].4s, %[sm1].4s, %[tp0].4s
    ext %[tab].16b, %[tab].16b, %[tab].16b, #4   // [b5, a5, a3, b3]
    zip1 %[tp0].2s, %[tab].2s, %[ta0].2s         // [b5, a0, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b5, b5, a0, a0]
    fmla %[dt2].4s, %[sm3].4s, %[tp0].4s
    fmls %[dt3].4s, %[sm2].4s, %[tp0].4s
    fdiv %[dt0].4s, %[dt0].4s, %[det].4s
    fdiv %[dt1].4s, %[dt1].4s, %[det].4s
    fdiv %[dt2].4s, %[dt2].4s, %[det].4s
    fdiv %[dt3].4s, %[dt3].4s, %[det].4s)"
    : [det]"=&x"(det), [tp0]"=&x"(t0), [tp1]"=&x"(t1)
    , [ta0]"=&x"(a0), [tb0]"=&x"(b0), [tab]"=&x"(ab)
    , [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    , [sm0]"=&x"(s0), [sm1]"=&x"(s1), [sm2]"=&x"(s2), [sm3]"=&x"(s3)
    : [mt0]"x"(c0), [mt1]"x"(c1), [mt2]"x"(c2), [mt3]"x"(c3));
    return {r0, r1, r2, r3};
  }


  inline bool inverse (f32 threshold = 0.0000001F) noexcept {
    f32v4 det, t0, t1, a0, b0, ab;
    f32 abs;
    asm volatile (R"(
    dup %[tp0].4s, %[mt0].s[0]                   // [s00, s00, s00, s00]
    ins %[tp0].s[3], %[mt0].s[1]                 // [s00, s00, s00, s01]
    ext %[tp1].16b, %[mt1].16b, %[mt1].16b, #4   // [s11, s12, s13, s10]
    ins %[tp1].s[3], %[mt1].s[3]                 // [s11, s12, s13, s13]
    fmul %[ta0].4s, %[tp0].4s, %[tp1].4s         // [a0, a1, a2, a4]
    ext %[tp0].16b, %[mt0].16b, %[mt0].16b       // [s01, s02, s03, s01]
    ins %[tp0].s[3], %[mt0].s[3]                 // [s01, s02, s03, s03]
    dup %[tp1].4s, %[mt1].s[0]                   // [s10, s10, s10, s10]
    ins %[tp1].s[3], %[mt1].s[1]                 // [s10, s10, s10, s11]
    fmls %[ta0].4s, %[tp0].4s, %[tp1].4s         // [a0, a1, a2, a4]
    dup %[tp0].4s, %[mt2].s[1]                   // [s21, s21, s21, s21]
    ins %[tp0].s[0], %[mt2].s[2]                 // [s22, s21, s21, s21]
    ins %[tp0].s[3], %[mt2].s[0]                 // [s22, s21, s21, s20]
    dup %[tp1].2s, %[mt3].s[3]                   // [s33, s33, -, -]
    ins %[tp1].s[2], %[mt3].s[2]                 // [s33, s33, s32, -]
    ins %[tp1].s[3], %[mt3].s[2]                 // [s33, s33, s32, s32]
    fmul %[tb0].4s, %[tp0].4s, %[tp1].4s         // [b0, b1, b2, b4]
    dup %[tp0].2s, %[mt2].s[3]                   // [s23, s23, -, -]
    ins %[tp0].s[2], %[mt2].s[2]                 // [s23, s23, s22, -]
    ins %[tp0].s[3], %[mt2].s[2]                 // [s23, s23, s22, s22]
    dup %[tp1].4s, %[mt3].s[1]                   // [s31, s31, s31, s31]
    ins %[tp1].s[0], %[mt3].s[2]                 // [s32, s31, s31, s31]
    ins %[tp1].s[3], %[mt3].s[0]                 // [s32, s31, s31, s30]
    fmls %[tb0].4s, %[tp0].4s, %[tp1].4s         // [b0, b1, b2, b4]
    ext %[tp0].16b, %[mt0].16b, %[mt2].16b, #4   // [s01, s02, s03, s20]
    ins %[tp0].s[2], %[mt2].s[0]                 // [s01, s02, s20, s20]
    ext %[tp1].16b, %[mt1].16b, %[mt3].16b, #8   // [s12, s13, s30, s31]
    ins %[tp1].s[2], %[mt3].s[3]                 // [s12, s13, s33, s31]
    fmul %[ta1].4s, %[tp0].4s, %[tp1].4s         // [a3, a5, b3, b5]
    ext %[tp0].16b, %[mt0].16b, %[mt2].16b, #8   // [s02, s03, s20, s21]
    ins %[tp0].s[2], %[mt2].s[3]                 // [s02, s03, s23, s21]
    ext %[tp1].16b, %[mt1].16b, %[mt3].16b, #4   // [s11, s12, s13, s30]
    ins %[tp1].s[2], %[mt3].s[0]                 // [s11, s12, s30, s30]
    fmls %[tab].4s, %[tp0].4s, %[tp1].4s         // [a3, a5, b3, b5]
    ins %[tp1].d[0], %[tab].d[1]                 // [b3, b5, -, -]
    fmul %[tp1].2s, %[tab].2s, %[tb1].2s         // [a3b3, a5b5, -, -]
    eor %[tp0].8b, %[tp0].8b, %[tp0].8b          // [0, 0, -, -]
    zip1 %[tp1].2d, %[tp1].2d, %[tp0].2d         // [a3b3, a5b5, 0, 0]
    fmul %[tp0].4s, %[ta0].4s, %[tb0].4s         // [a0b0, a1b1, a2b2, a4b4]
    uzp1 %[det].4s, %[tp0].4s, %[tp1].4s         // [a0b0, a2b2, a3b3, 0]
    uzp2 %[tp1].4s, %[tp0].4s, %[tp1].4s         // [a1b1, a4b4, a5b5, 0]
    fneg %[tp1].2s, %[tp1].2s                    // [-a1b1, -a4b4, a5b5, 0]
    fadd %[det].4s, %[det].4s, %[tp1].4s         // [a0b0-a1b1, a2b2-a4b4, a3b3+a5b5, 0]
    faddp %[det].4s, %[det].4s, %[det].4s
    faddp %[det].4s, %[det].4s, %[det].4s
    fabs %s[abs], %s[det])"
    : [abs]"=x"(abs), [det]"=x"(det), [tp0]"=&x"(t0), [tp1]"=&x"(t1)
    , [ta0]"=&x"(a0), [tb0]"=&x"(b0), [tab]"=&x"(ab)
    : [mt0]"x"(c0), [mt1]"x"(c1), [mt2]"x"(c2), [mt3]"x"(c3));
    if (abs < threshold)
      return false;
    f32v4 s0, s1, s2, s3;
    asm volatile (R"(
    fneg %[mt0].4s, %[mt0].4s                    // [-s00, -s01, -s02, -s03]
    fneg %[mt2].4s, %[mt2].4s                    // [-s20, -s21, -s22, -s23]
    zip1 %[tp0].2d, %[mt0].2d, %[mt2].2d         // [-s00, -s01, -s20, -s21]
    zip1 %[tp1].2d, %[mt1].2d, %[mt3].2d         // [s10, s11, s30, s31]
    trn1 %[sm0].4s, %[tp1].4s, %[tp0].4s         // [s10, -s00, s30, -s20]
    trn2 %[sm1].4s, %[tp1].4s, %[tp0].4s         // [s11, -s01, s31, -s21]
    zip2 %[tp0].2d, %[mt0].2d, %[mt2].2d         // [-s02, -s03, -s22, -s23]
    zip2 %[tp1].2d, %[mt1].2d, %[mt3].2d         // [s12, s13, s32, s33]
    trn1 %[sm2].4s, %[tp1].4s, %[tp0].4s         // [s12, -s02, s32, -s22]
    trn2 %[sm3].4s, %[tp1].4s, %[tp0].4s         // [s13, -s03, s33, -s23]
    eor %[mt1].16b, %[mt1].16b, %[mt1].16b       // [0, 0, 0, 0]
    mov %[mt3].16b, %[mt1].16b                   // [0, 0, 0, 0]
    rev64 %[tab].2s, %[tab].2s                   // [a5, a3, b3, b5]
    zip1 %[tp0].2s, %[tb0].2s, %[tab].2s         // [b0, a5, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b0, b0, a5, a5]
    fmul %[mt0].4s, %[sm1].4s, %[tp0].4s
    fmls %[mt1].4s, %[sm0].4s, %[tp0].4s
    ext %[tp1].16b, %[ta0].16b, %[ta0].16b, #8   // [a2, a4, a0, a1]
    zip2 %[tp0].2s, %[tb0].2s, %[tp1].2s         // [b1, a4, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b1, b1, a4, a4]
    fmls %[mt0].4s, %[sm2].4s, %[tp0].4s
    fmul %[mt2].4s, %[sm0].4s, %[tp0].4s
    ext %[tb0].16b, %[tb0].16b, %[tb0].16b, #4   // [b1, b2, b4, b0]
    zip2 %[tp0].2s, %[tb0].2s, %[tab].2s         // [b2, a3, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b2, b2, a3, a3]
    fmla %[mt0].4s, %[sm3].4s, %[tp0].4s
    fmls %[mt3].4s, %[sm0].4s, %[tp0].4s
    ext %[tab].16b, %[tab].16b, %[tab].16b, #8   // [b3, b5, a5, a3]
    zip1 %[tp0].2s, %[tab].2s, %[tp1].2s         // [b3, a2, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b3, b3, a2, a2]
    fmla %[mt1].4s, %[sm2].4s, %[tp0].4s
    fmls %[mt2].4s, %[sm1].4s, %[tp0].4s
    ext %[tb0].16b, %[tb0].16b, %[tb0].16b, #4   // [b2, b4, b0, b1]
    zip2 %[tp0].2s, %[tb0].2s, %[ta0].2s         // [b4, a1, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b4, b4, a1, a1]
    fmls %[mt1].4s, %[sm3].4s, %[tp0].4s
    fmla %[mt3].4s, %[sm1].4s, %[tp0].4s
    ext %[tab].16b, %[tab].16b, %[tab].16b, #4   // [b5, a5, a3, b3]
    zip1 %[tp0].2s, %[tab].2s, %[ta0].2s         // [b5, a0, -, -]
    zip1 %[tp0].4s, %[tp0].4s, %[tp0].4s         // [b5, b5, a0, a0]
    fmla %[mt2].4s, %[sm3].4s, %[tp0].4s
    fmls %[mt3].4s, %[sm2].4s, %[tp0].4s
    fdiv %[mt0].4s, %[mt0].4s, %[det].4s
    fdiv %[mt1].4s, %[mt1].4s, %[det].4s
    fdiv %[mt2].4s, %[mt2].4s, %[det].4s
    fdiv %[mt3].4s, %[mt3].4s, %[det].4s)"
    : [tp0]"=&x"(t0), [tp1]"=&x"(t1)
    , [sm0]"=&x"(s0), [sm1]"=&x"(s1), [sm2]"=&x"(s2), [sm3]"=&x"(s3)
    , [ta0]"+x"(a0), [tb0]"+x"(b0), [tab]"+x"(ab)
    , [mt0]"+x"(c0), [mt1]"+x"(c1), [mt2]"+x"(c2), [mt3]"+x"(c3)
    : [det]"x"(det));
    return true;
  }


  inline static Mat4 zero () noexcept {
    f32v4 r;
    asm volatile (R"(
    eor %[dst].16b, %[dst].16b, %[dst].16b)"
    : [dst]"=r"(r));
    return {r, r, r, r};
  }
  inline static Mat4 identity () noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    mov %[dt3].16b, %[dt0].16b
    ins %[dt0].s[0], %[imm].s[0]
    ins %[dt1].s[1], %[imm].s[0]
    ins %[dt2].s[2], %[imm].s[0]
    ins %[dt3].s[3], %[imm].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [imm]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 translate (Vec3 p) noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    ins %[dt0].s[0], %[one].s[0]
    ins %[dt1].s[1], %[one].s[0]
    ins %[dt2].s[2], %[one].s[0]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=x"(r3)
    : "[dt3]"(p.p), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 rotatex (f32 a) noexcept {
    auto s = sinf(a);
    auto c = cosf(a);
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    ins %[dt1].s[1], %[cos].s[0]
    ins %[dt2].s[2], %[cos].s[0]
    ins %[dt1].s[2], %[sin].s[0]
    fneg %s[dt3], %s[sin]
    ins %[dt2].s[1], %[dt3].s[0]
    mov %[dt3].16b, %[dt0].16b
    ins %[dt0].s[0], %[one].s[0]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [sin]"x"(s), [cos]"x"(c), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 rotatey (f32 a) noexcept {
    auto s = sinf(a);
    auto c = cosf(a);
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    ins %[dt0].s[0], %[cos].s[0]
    ins %[dt2].s[2], %[cos].s[0]
    ins %[dt2].s[0], %[sin].s[0]
    fneg %s[dt3], %s[sin]
    ins %[dt0].s[2], %[dt3].s[0]
    mov %[dt3].16b, %[dt1].16b
    ins %[dt1].s[1], %[one].s[0]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [sin]"x"(s), [cos]"x"(c), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 rotatez (f32 a) noexcept {
    auto s = sinf(a);
    auto c = cosf(a);
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    ins %[dt0].s[0], %[cos].s[0]
    ins %[dt1].s[1], %[cos].s[0]
    ins %[dt0].s[1], %[sin].s[0]
    fneg %s[dt3], %s[sin]
    ins %[dt1].s[0], %[dt3].s[0]
    mov %[dt3].16b, %[dt2].16b
    ins %[dt2].s[2], %[one].s[0]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [sin]"x"(s), [cos]"x"(c), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 rotatexp (f32 a, Vec3 p) noexcept {
    auto s = sinf(a);
    auto c = cosf(a);
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    ins %[dt1].s[1], %[cos].s[0]
    ins %[dt2].s[2], %[cos].s[0]
    ins %[dt1].s[2], %[sin].s[0]
    fneg %s[dt3], %s[sin]
    ins %[dt2].s[1], %[dt3].s[0]
    mov %[dt3].16b, %[dot].16b
    fmls %[dt3].4s, %[dt1].4s, %[dot].s[1]
    fmls %[dt3].4s, %[dt2].4s, %[dot].s[2]
    ins %[dt3].s[0], %[dt1].s[0]
    ins %[dt0].s[0], %[one].s[0]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [dot]"x"(p.p), [sin]"x"(s), [cos]"x"(c), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 rotateyp (f32 a, Vec3 p) noexcept {
    auto s = sinf(a);
    auto c = cosf(a);
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    ins %[dt0].s[0], %[cos].s[0]
    ins %[dt2].s[2], %[cos].s[0]
    ins %[dt2].s[0], %[sin].s[0]
    fneg %s[dt3], %s[sin]
    ins %[dt0].s[2], %[dt3].s[0]
    mov %[dt3].16b, %[dot].16b
    fmls %[dt3].4s, %[dt0].4s, %[dot].s[0]
    fmls %[dt3].4s, %[dt2].4s, %[dot].s[2]
    ins %[dt3].s[1], %[dt0].s[1]
    ins %[dt1].s[1], %[one].s[0]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [dot]"x"(p.p), [sin]"x"(s), [cos]"x"(c), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 rotatezp (f32 a, Vec3 p) noexcept {
    auto s = sinf(a);
    auto c = cosf(a);
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    ins %[dt0].s[0], %[cos].s[0]
    ins %[dt1].s[1], %[cos].s[0]
    ins %[dt0].s[1], %[sin].s[0]
    fneg %s[dt3], %s[sin]
    ins %[dt1].s[0], %[dt3].s[0]
    mov %[dt3].16b, %[dot].16b
    fmls %[dt3].4s, %[dt0].4s, %[dot].s[0]
    fmls %[dt3].4s, %[dt1].4s, %[dot].s[1]
    ins %[dt3].s[2], %[dt0].s[2]
    ins %[dt2].s[2], %[one].s[0]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [dot]"x"(p.p), [sin]"x"(s), [cos]"x"(c), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 rotatev (f32 a, Vec3 v) noexcept {
    auto s = sinf(a);
    auto c = 1.0F - cosf(a);
    f32v4 r0, r1, r2, r3, t;
    asm volatile (R"(
    ins %[dt0].s[0], %[cos].s[0]
    ins %[dt1].s[1], %[cos].s[1]
    ins %[dt2].s[2], %[cos].s[2]
    fmul %[dt3].4s, %[vec].4s, %[sin].s[0]
    fneg %[tmp].4s, %[dt3].4s
    ins %[dt0].s[1], %[dt3].s[2]
    ins %[dt0].s[2], %[tmp].s[1]
    ins %[dt1].s[0], %[tmp].s[2]
    ins %[dt1].s[2], %[dt3].s[0]
    ins %[dt2].s[0], %[dt3].s[1]
    ins %[dt2].s[1], %[tmp].s[0]
    fmul %[tmp].4s, %[vec].4s, %[cos].s[0]
    fmla %[dt0].4s, %[tmp].4s, %[vec].s[0]
    fmla %[dt1].4s, %[tmp].4s, %[vec].s[1]
    fmla %[dt2].4s, %[tmp].4s, %[vec].s[2]
    eor %[dt3].16b, %[dt3].16b, %[dt3].16b
    ins %[dt0].s[3], %[dt3].s[0]
    ins %[dt1].s[3], %[dt3].s[0]
    ins %[dt2].s[3], %[dt3].s[0]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3), [tmp]"=&x"(t)
    : [vec]"x"(v.p), [sin]"x"(s), [cos]"x"(c), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 rotatevp (f32 a, Vec3 v, Vec3 p) noexcept {
    auto s = sinf(a);
    auto c = 1.0F - cosf(a);
    f32v4 r0, r1, r2, r3, t;
    asm volatile (R"(
    ins %[dt0].s[0], %[cos].s[0]
    ins %[dt1].s[1], %[cos].s[1]
    ins %[dt2].s[2], %[cos].s[2]
    fmul %[dt3].4s, %[vec].4s, %[sin].s[0]
    fneg %[tmp].4s, %[dt3].4s
    ins %[dt0].s[1], %[dt3].s[2]
    ins %[dt0].s[2], %[tmp].s[1]
    ins %[dt1].s[0], %[tmp].s[2]
    ins %[dt1].s[2], %[dt3].s[0]
    ins %[dt2].s[0], %[dt3].s[1]
    ins %[dt2].s[1], %[tmp].s[0]
    fmul %[tmp].4s, %[vec].4s, %[cos].s[0]
    fmla %[dt0].4s, %[tmp].4s, %[vec].s[0]
    fmla %[dt1].4s, %[tmp].4s, %[vec].s[1]
    fmla %[dt2].4s, %[tmp].4s, %[vec].s[2]
    eor %[dt3].16b, %[dt3].16b, %[dt3].16b
    ins %[dt0].s[3], %[dt3].s[0]
    ins %[dt1].s[3], %[dt3].s[0]
    ins %[dt2].s[3], %[dt3].s[0]
    mov %[dt3].16b, %[dot].16b
    fmls %[dt3].4s, %[dt0].4s, %[dot].s[0]
    fmls %[dt3].4s, %[dt1].4s, %[dot].s[1]
    fmls %[dt3].4s, %[dt2].4s, %[dot].s[2]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3), [tmp]"=&x"(t)
    : [vec]"x"(v.p), [dot]"x"(p.p), [sin]"x"(s), [cos]"x"(c), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 scale (Vec3 t) noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    mov %[dt3].16b, %[dt0].16b
    ins %[dt0].s[0], %[arg].s[0]
    ins %[dt1].s[1], %[arg].s[1]
    ins %[dt2].s[2], %[arg].s[2]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [arg]"x"(t.p), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 scalep (Vec3 t, Vec3 p) noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    mov %[dt3].16b, %[dot].16b
    ins %[dt0].s[0], %[arg].s[0]
    ins %[dt1].s[1], %[arg].s[1]
    ins %[dt2].s[2], %[arg].s[2]
    fmls %[dt3].4s, %[arg].4s, %[dot].4s
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [arg]"x"(t.p), [dot]"x"(p.p), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 scalev (f32 t, Vec3 v) noexcept {
    f32v4 r0, r1, r2, r3, t;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    mov %[dt3].16b, %[dt0].16b
    ins %[dt0].s[0], %[arg].s[0]
    ins %[dt1].s[1], %[arg].s[1]
    ins %[dt2].s[2], %[arg].s[2]
    fmul %[tmp].4s, %[vec].4s, %[arg].s[0]
    ins %[tmp].s[3], %[dt3].s[0]
    fmla %[dt0].4s, %[tmp].4s, %[vec].s[0]
    fmla %[dt1].4s, %[tmp].4s, %[vec].s[1]
    fmla %[dt2].4s, %[tmp].4s, %[vec].s[2]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3), [tmp]"=&x"(t)
    : [arg]"x"(t-1.0F), [vec]"x"(v.p), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 scalevp (f32 t, Vec3 v, Vec3 p) noexcept {
    f32v4 r0, r1, r2, r3, t;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    mov %[dt3].16b, %[dt0].16b
    fmul %[tp0].4s, %[vec].4s, %[dot].4s
    faddp %[tp0].2s, %[tp0].2s, %[tp0].2s
    ins %[tp1].s[0], %[tp0].s[2]
    fadd %s[tp0], %s[tp0], %s[tp1]
    fmul %[tp1].4s, %[vec].4s, %[arg].s[0]
    ins %[tp1].s[3], %[dt3].s[0]
    ins %[dt0].s[0], %[arg].s[0]
    ins %[dt1].s[1], %[arg].s[1]
    ins %[dt2].s[2], %[arg].s[2]
    fmla %[dt0].4s, %[tmp].4s, %[vec].s[0]
    fmla %[dt1].4s, %[tmp].4s, %[vec].s[1]
    fmla %[dt2].4s, %[tmp].4s, %[vec].s[2]
    fmls %[dt3].4s, %[tmp].4s, %[tp0].s[0]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3), [tp0]"=&x"(t0), [tp1]"=&x"(t1)
    : [arg]"x"(t-1.0F), [vec]"x"(v.p), [dot]"x"(p.p), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 skewx (Vec3 v) noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    dup %[dt0].4s, %[vec].s[0]
    fdiv %[dt0].4s, %[vec].4s, %[dt0].4s
    eor %[dt3].16b, %[dt3].16b, %[dt3].16b
    ins %[dt0].s[3], %[dt3].s[0]
    mov %[dt1].16b, %[dt3].16b
    mov %[dt2].16b, %[dt3].16b
    ins %[dt1].s[1], %[dt0].s[0]
    ins %[dt2].s[2], %[dt0].s[0]
    ins %[dt3].s[3], %[dt0].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [vec]"x"(v.p));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 skewy (Vec3 v) noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    dup %[dt1].4s, %[vec].s[1]
    fdiv %[dt1].4s, %[vec].4s, %[dt1].4s
    eor %[dt3].16b, %[dt3].16b, %[dt3].16b
    ins %[dt1].s[3], %[dt3].s[0]
    mov %[dt0].16b, %[dt3].16b
    mov %[dt2].16b, %[dt3].16b
    ins %[dt0].s[0], %[dt1].s[1]
    ins %[dt2].s[2], %[dt1].s[1]
    ins %[dt3].s[3], %[dt1].s[1])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [vec]"x"(v.p));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 skewz (Vec3 v) noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    dup %[dt2].4s, %[vec].s[2]
    fdiv %[dt2].4s, %[vec].4s, %[dt2].4s
    eor %[dt3].16b, %[dt3].16b, %[dt3].16b
    ins %[dt2].s[3], %[dt3].s[0]
    mov %[dt0].16b, %[dt3].16b
    mov %[dt1].16b, %[dt3].16b
    ins %[dt0].s[0], %[dt2].s[2]
    ins %[dt1].s[1], %[dt2].s[2]
    ins %[dt3].s[3], %[dt2].s[2])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [vec]"x"(v.p));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 skewxp (Vec3 v, f32 p) noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    dup %[dt0].4s, %[vec].s[0]
    fdiv %[dt0].4s, %[vec].4s, %[dt0].4s
    eor %[dt3].16b, %[dt3].16b, %[dt3].16b
    ins %[dt0].s[3], %[dt3].s[0]
    mov %[dt1].16b, %[dt3].16b
    mov %[dt2].16b, %[dt3].16b
    ins %[dt1].s[1], %[dt0].s[0]
    ins %[dt2].s[2], %[dt0].s[0]
    fmls %[dt3].4s, %[dt0].4s, %[dot].4s
    ins %[dt3].s[0], %[dt0].s[3]
    ins %[dt3].s[3], %[dt0].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [vec]"x"(v.p), [dot]"x"(p.p));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 skewyp (Vec3 v, f32 p) noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    dup %[dt1].4s, %[vec].s[1]
    fdiv %[dt1].4s, %[vec].4s, %[dt1].4s
    eor %[dt3].16b, %[dt3].16b, %[dt3].16b
    ins %[dt1].s[3], %[dt3].s[0]
    mov %[dt0].16b, %[dt3].16b
    mov %[dt2].16b, %[dt3].16b
    ins %[dt0].s[0], %[dt1].s[1]
    ins %[dt2].s[2], %[dt1].s[1]
    fmls %[dt3].4s, %[dt1].4s, %[dot].4s
    ins %[dt3].s[1], %[dt1].s[3]
    ins %[dt3].s[3], %[dt1].s[1])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [vec]"x"(v.p), [dot]"x"(p.p));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 skewzp (Vec3 v, f32 p) noexcept {
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    dup %[dt2].4s, %[vec].s[2]
    fdiv %[dt2].4s, %[vec].4s, %[dt2].4s
    eor %[dt3].16b, %[dt3].16b, %[dt3].16b
    ins %[dt2].s[3], %[dt3].s[0]
    mov %[dt0].16b, %[dt3].16b
    mov %[dt1].16b, %[dt3].16b
    ins %[dt0].s[0], %[dt2].s[2]
    ins %[dt1].s[1], %[dt2].s[2]
    fmls %[dt3].2s, %[dt2].2s, %[dot].2s
    ins %[dt3].s[3], %[dt2].s[2])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [vec]"x"(v.p), [dot]"x"(p.p));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 lookat (Vec3 s, Vec3 d, Vec3 u) noexcept {
    Vec3 z = s - d;
    z.normalize();
    Vec3 x = u ^ z;
    x.normalize();
    Vec3 y = z ^ x;
    y.normalize();
    f32v4 r0, r1, r2, r3, t0, t1;
    asm volatile (R"(
    eor %[dt3].16b, %[dt3].16b, %[dt3].16b
    zip1 %[tp0].2d, %[rvx].2d, %[rvz].2d
    zip1 %[tp1].2d, %[rvy].2d, %[dt3].2d
    trn1 %[dt0].4s, %[tp0].4s, %[tp1].4s
    trn2 %[dt1].4s, %[tp0].4s, %[tp1].4s
    zip2 %[tp0].2d, %[rvx].2d, %[rvz].2d
    zip2 %[tp1].2d, %[rvy].2d, %[dt3].2d
    trn1 %[dt2].4s, %[tp0].4s, %[tp1].4s
    eor %[dt3].16b, %[dt3].16b, %[dt3].16b
    fmls %[dt3].4s, %[dt0].4s, %[src].s[0]
    fmls %[dt3].4s, %[dt1].4s, %[src].s[1]
    fmls %[dt3].4s, %[dt2].4s, %[src].s[2]
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3), [tp0]"=&x"(t0), [tp1]"=&x"(t1)
    : [rvx]"x"(x.p), [rvy]"x"(y.p), [rvz]"x"(z.p), [src]"x"(s.p), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 orthogonal (const Box3& frustum) noexcept {
    auto dv = frustum.max - frustum.min;
    auto av = frustum.min + frustum.max;
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    fadd %s[dt3], %s[one], %s[one]
    dup %[dt3].4s, %[dt3].s[0]
    fdiv %[dt3].4s, %[dt3].4s, %[dff].4s
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    ins %[dt0].s[0], %[dt3].s[0]
    ins %[dt1].s[1], %[dt3].s[1]
    ins %[dt2].s[2], %[dt3].s[2]
    fdiv %[dt3].4s, %[sum].4s, %[dff].4s
    ins %[dt3].s[3], %[one].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [sum]"x"(av), [dff]"x"(dv), [one]"x"(1.0F));
    return {r0, r1, r2, r3};
  }
  inline static Mat4 perspective (f32 fovy, f32 aspect, f32 near, f32 far) noexcept {
    auto i = 1.0F / (near-far);
    auto a = 1.0F / tanf(0.5F*fovy);
    auto b = a / aspect;
    auto c = i * (far+near);
    auto d = -2.0F * i * far * near;
    f32v4 r0, r1, r2, r3;
    asm volatile (R"(
    eor %[dt0].16b, %[dt0].16b, %[dt0].16b
    mov %[dt1].16b, %[dt0].16b
    mov %[dt2].16b, %[dt0].16b
    mov %[dt3].16b, %[dt0].16b
    ins %[dt0].s[0], %[ppb].s[0]
    ins %[dt1].s[1], %[ppa].s[0]
    ins %[dt2].s[2], %[ppc].s[0]
    ins %[dt2].s[3], %[one].s[0]
    ins %[dt3].s[2], %[ppd].s[0])"
    : [dt0]"=&x"(r0), [dt1]"=&x"(r1), [dt2]"=&x"(r2), [dt3]"=&x"(r3)
    : [ppa]"x"(a), [ppb]"x"(b), [ppc]"x"(c), [ppd]"x"(d), [one]"x"(-1.0F));
    return {r0, r1, r2, r3};
  }
};


template<>
inline f32 Mat4::c<0> () const noexcept {
  return c0[0];
}
template<>
inline f32 Mat4::c<1> () const noexcept {
  return c0[1];
}
template<>
inline f32 Mat4::c<2> () const noexcept {
  return c0[2];
}
template<>
inline f32 Mat4::c<3> () const noexcept {
  return c0[3];
}
template<>
inline f32 Mat4::c<4> () const noexcept {
  return c1[0];
}
template<>
inline f32 Mat4::c<5> () const noexcept {
  return c1[1];
}
template<>
inline f32 Mat4::c<6> () const noexcept {
  return c1[2];
}
template<>
inline f32 Mat4::c<7> () const noexcept {
  return c1[3];
}
template<>
inline f32 Mat4::c<8> () const noexcept {
  return c2[0];
}
template<>
inline f32 Mat4::c<9> () const noexcept {
  return c2[1];
}
template<>
inline f32 Mat4::c<10> () const noexcept {
  return c2[2];
}
template<>
inline f32 Mat4::c<11> () const noexcept {
  return c2[3];
}
template<>
inline f32 Mat4::c<12> () const noexcept {
  return c3[0];
}
template<>
inline f32 Mat4::c<13> () const noexcept {
  return c3[1];
}
template<>
inline f32 Mat4::c<14> () const noexcept {
  return c3[2];
}
template<>
inline f32 Mat4::c<15> () const noexcept {
  return c3[3];
}


inline Vec4 operator* (Vec4 v, const Mat4& m) noexcept {
  return ~m * v;
}


CLOSE_JLIB_NS